#include <iostream>
#include <vector>
#include <ctime>
#include <numeric>

#include "tensor.h"
#include "covariance.h"
#include "mtrand.h"
#include "dataset.h"
#include "pdf.h"
#include <vector>
#include <sstream>
#include <fstream>

#include "TFile.h"
#include "TChain.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TProfile2D.h"
#include "TF1.h"
#include "TF2.h"
#include "TF3.h"
#include "TLorentzVector.h"
#include "TGraphErrors.h"
#include "TGraph2DErrors.h"
#include "make_list.h"

using namespace std;

//initialize random number generator
MTRand_int32 irand(time(0));

const double PI=3.14159265359;
const double min_pt=320;
const double tau32_cut=0.4;
const double mass_cut=300;

//fast_mode skips template computatations to reduce memory usage
const bool fast_mode=true;

//counter for number of events
class Counter{
  
public:
  double n;
  double n2;
  double error;
  double n_corr;
  double n2_corr;

  void Add(double weight, double error_=0)
  { 
    n+=weight;
    n2+=weight*weight;
    
    error+=error_;

    double corr_weight=weight-error_;
    n_corr+=corr_weight;
    n2_corr+=corr_weight*corr_weight;
  }

  double nevt() const
  {
    return n;
  }

  double nevt_corr() const
  {
    return n_corr;
  }
    
  double err() const
  {
    return sqrt(n2 + error*error);
  }

  double bias() const
  {
    return error;
  }

  double var() const
  {
    return sqrt(n2);
  }

  double var_corr() const
  {
    return sqrt(n2_corr);
  }

};

class ABCD_Counter
{
public:
  double t32_cut;
  double mj_cut;

  double a,b,c,d;
  double a2,b2,c2,d2;

  ABCD_Counter(double t32, double mj):
    t32_cut(t32), mj_cut(mj)
  {}
  
  void Add(double mj1, double mj2, 
	   double t32_1, double t32_2, double w)
  {
    double w2=w*w;

    if(mj1<mj_cut && mj2<mj_cut)
      {
	a+=w;
	a2+=w2;
      }

    else if(mj1>mj_cut && t32_1<t32_cut &&
	    mj2<mj_cut)
      {
	b+=w;
	b2+=w2;
      }

    else if(mj2>mj_cut && t32_2<t32_cut &&
	    mj1<mj_cut)
      {
	c+=w;
	c2+=w2;
      }

    else if(mj1>mj_cut && t32_cut>t32_1 &&
	    mj2>mj_cut && t32_cut>t32_2 )
      {
	d+=w;
	d2+=w2;
      }
  }

  
  double factor() const
  {
    if(b==0 || c==0) return 0;
    
    return d*a/(b*c);
  }

  double factor_err() const
  {
    double num[4]={a,b,c,d};
    double num2[4]={a2,b2,c2,d2};
    
    double result_err=0;
    for(unsigned int i=0; i<4; ++i)
      {
	double temp_err=0;
	if(num[i]>0)
	  temp_err=num2[i]/(num[i]*num[i]);

	result_err+=temp_err;
      }

    //return full error
    return factor()*result_err;
  }






};




//function to produce a 2D plot from DataSet
TGraph2D plot2D(DataSet& data, char* name)
{
  //ensure that dataset is two dimensional
  assert(data.Rank() == 2);

  PDF pdf=data.ComputePDF();
  
  //declare some x,y points
  vector<double> x;
  vector<double> y;
  vector<double> val;
  
  const int nbin=20;

  double x_diff=(data.max[0]-data.min[0])/nbin;
  double y_diff=(data.max[1]-data.min[1])/nbin;

  for(unsigned int i=0; i<nbin; ++i)
  for(unsigned int j=0; j<nbin; ++j)
    {
      double x_val=data.min[0]+x_diff*(i+0.5);
      double y_val=data.min[1]+y_diff*(j+0.5);

      x.push_back(x_val);
      y.push_back(y_val);
      
      val.push_back(pdf(make_list<double>(x_val, y_val)));
    }
    
  //get a new TGraph2D
  TGraph2D graph(x.size(), &x[0], &y[0], &val[0]);
  graph.SetName(name);
  graph.SetTitle(name);

  return graph;
}


class Vec
{
public:
  vector<double> x,y,err;

  Vec(){}

  void push_back(double x_, double y_, double yerr_)
  {
    x.push_back(x_);
    y.push_back(y_);
    err.push_back(yerr_);    
  }

  TGraphErrors get_graph(string name, bool logx=false) const
  {
    vector<double> xprime(x);
    if(logx)
      {
	for(unsigned int i=0; i<xprime.size(); ++i)
	  {
	    xprime[i]=exp10(-2*xprime[i]);
	  }
      }
    TGraphErrors result(xprime.size(), &xprime[0], &y[0], NULL, &err[0]);
    result.SetTitle(name.c_str());
    result.SetName(name.c_str());

    return result;
  }

  //get a 1D graph from pt
  void fill_data(PDF& pdf, double pt)
  {
    PDF::MC_Generator gen=
      pdf.Generator(make_list<PDF::MC_Flag>
		    (PDF::OUTPUT, PDF::INPUT));

    vector<PDF::MC_Event> evt=
      gen.GenerateFull(make_list<double>(log(pt/min_pt)));

    for(unsigned int i=0; i<evt.size(); ++i)
      {
	x.push_back(evt[i].value[0]);
	y.push_back(evt[i].weight);
	err.push_back(evt[i].err);
      }
  }

};



float DeltaR(float eta, float phi)
{
  if(fabs(phi)> PI)
    phi= PI-fabs(phi);
  return sqrt(eta*eta + phi*phi);
}

int main(int argc, char **argv)
{

  //data-driven analysis on RPV-gluinos 1210.4813
  if(argc < 2+1)
    {
      cout<<"Usage: ./RPV_gluino madgraph_input output.root [extra smoothing] [random seed if PDF_smearing on] <LHC> <LHC lumi in pb^{-1}>"
	  <<endl;
      throw;
    }

  cout<<"madgraph inputs are read as madgraph_input.?.root ranges from 1 to 12"
      <<endl;
  cout<<"Outputs are automatically stored in PDF_*"<<endl;

  //see if we want to fluctuate the PDF to get variance
  bool b_smear_PDF=false;
  bool b_smear_lhc=false;
  double lhc_factor=100.0*1e9;


  //boolean to divide the job
  bool b_divide=false;
  int n_div=100;
  int n_part=0;

  if(argc >= 3+1)
    {
      cout<<"Extra smoothing turned on: "<<atof(argv[3])<<endl;
      DataSet::kernel_smoothing_factor=atof(argv[3]);
    }

  if(argc >= 4+1)
    {
      b_smear_PDF=true;

      cout<<"PDF smearing on, with seed: "<<atoi(argv[4])<<endl;

      //set random seed
      DataSet::SetRandomSeed(atoi(argv[4]));


      //see if there are extra arguments
      if(argc >= 6+1 && string(argv[5])=="LHC")
	{
	  b_smear_lhc=true;
	  cout<<"LHC mode enabled"<<endl;

	  //lumi in picobarn
	  lhc_factor= atof(argv[6])*1e9;

	  cout<<"scale up to "<<atof(argv[6])<<" inverse picobarn"<<endl;


	  //which part to process
	  if(argc >= 7+1)
	    {
	      b_divide=true;
	      n_part=atoi(argv[7]);
	      cout<<"N_division enabled: "<<n_div<<" part "
		  <<n_part<<endl;	      

	      assert(n_part<n_div && n_part >=0);
	    }
	}
    }


  //needs to be changed
  //b_smear_PDF=false;
  //b_smear_lhc=false;

  //cout<<"WARNGING WARNING: QUICK CODE CHANGE LINE 329, NEEDS CHANGING"<<endl;
  //cout<<"WARNGING WARNING: OVERRIDE SETTINGS, NO SMEARING AND LHC"<<endl;


  //declare variables
  const int maxsize=20;
  Int_t           size;
  Double_t         w;   //[size]
  Float_t         pt[maxsize];   //[size]
  Float_t         eta[maxsize];
  Float_t         phi[maxsize];
  Float_t         m[maxsize];   //[size]
  Float_t         t1[maxsize];
  Float_t         t2[maxsize];
  Float_t         t3[maxsize];

  //list of DataSet
  vector<DataSet*> dataset_mpt_njet(1);
  vector<DataSet*> dataset_tau32_njet(1);
  vector<DataSet*> dataset_mpt_tau32_njet(1);

  //Statistics plots for different Pt
  TH2D stat_template_1("stat_template_1", "Template Stats (50, 136)",
		       100, 0, 3, 100, 0, 1.0);

  TH2D stat_template_2("stat_template_2", "Template Stats (136, 370)",
		       100, 0, 3, 100, 0, 1.0);

  TH2D stat_template_3("stat_template_3", "Template Stats (370, 1004)",
		       100, 0, 3, 100, 0, 1.0);

  TH2D stat_template_4("stat_template_4", "Template Stats (1004, 2730)",
		       100, 0, 3, 100, 0, 1.0);

  TH2D stat_template_5("stat_template_5", "Template Stats (2730, 7421)",
		       100, 0, 3, 100, 0, 1.0);

  
  TH2D* lhc_stat[5]={&stat_template_1, &stat_template_2, &stat_template_3,
		     &stat_template_4, &stat_template_5};

  vector<Counter> lhc_prescale(3, Counter());

  //no need to compute prescale
  double lhc_prescale_factor[3]={1.0, 1.0, 1.0};
  double lhc_prescale_cut[3]={100., 160., 250.};
  const double lhc_Hz=200.;


  //scale up to 20 fb^{-1}
  //weights are in milibarns, hence scale up by 1e12 to femtobarn
  int lhc_size=sizeof(lhc_stat)/sizeof(TH2D*);
  
  for(int is=0; is<lhc_size; ++is)
    {
      lhc_stat[is]->Sumw2();
    }
  
  //only do templates for two jets
  for(int n=0; n<2; n++)
    {
      if(!fast_mode)
	{
	  dataset_mpt_njet[n]=new
	    DataSet(make_list<double>(0, 0),
		    make_list<double>(2.5, 2.0),
		    make_list<double>(0.01, 0.01));
	  
	  dataset_tau32_njet[n]=new
	    DataSet(make_list<double>(0, 0),
		    make_list<double>(1.0, 2.0),
		    make_list<double>(0.01, 0.01));
	}

      if(n==1)
	dataset_mpt_tau32_njet[n]=dataset_mpt_tau32_njet[n-1];

      else
	dataset_mpt_tau32_njet[n]=new
	  DataSet(make_list<double>(0, 0, 0),
		  make_list<double>(2.5, 1.0, 2.0),
		  make_list<double>(0.01, 0.005, 0.008));
    }

  //smoothed plot
  DataSet* data_m12;

  if(!fast_mode)
  data_m12=new
    DataSet(make_list<double>(200, 200, 0, 0),
	    make_list<double>(2200, 2200, 1.0, 1.0),
	    make_list<double>(20, 20, 0.05, 0.05));


  //declare counters for real events and MC events
  vector<double> t32_cut;
  vector<double> mj_cut;
  vector<Counter> nevt_real;
  vector<ABCD_Counter> nevt_ABCD;
  vector<Counter> nevt_MC;

  //now do t32 cut from 0.2 to 1.0 in 0.1 increments (<- no cut)
  //jet mass cut from 0 to 500 GeV in 50 GeV increments
  for(double t32cut=0.2; t32cut<1.01; t32cut+=0.1)
    for(double mjcut=120; mjcut<1000; mjcut+=10.)
      {
	t32_cut.push_back(t32cut);
	mj_cut.push_back(mjcut);
	nevt_real.push_back(Counter());
	nevt_ABCD.push_back(ABCD_Counter(t32cut, mjcut));	
	nevt_MC.push_back(Counter());	
      }
  

  //now process madgraph events
  char alist[3]={'A', 'B', 'C'};



  //first get the prescale factors
  //do not do prescale

  /*
  if(false)
    {

  for(unsigned int k=1; k<=12; k++)
    for(unsigned int ka=0; ka<3; ka++)
      {
	TChain chain_madgraph("Ntuple");
	stringstream sout;
	sout<<argv[1]<<alist[ka]<<"."<<k<<".root";
	chain_madgraph.Add(sout.str().c_str());
	
	//read the input file
	chain_madgraph.SetBranchAddress("size", &size);
	chain_madgraph.SetBranchAddress("w", &w);
	chain_madgraph.SetBranchAddress("pt", pt);
	chain_madgraph.SetBranchAddress("eta", eta);
	chain_madgraph.SetBranchAddress("phi", phi);
	chain_madgraph.SetBranchAddress("m", m);
	chain_madgraph.SetBranchAddress("tau1", t1);
	chain_madgraph.SetBranchAddress("tau2", t2);
	chain_madgraph.SetBranchAddress("tau3", t3);
	
	double rescale=1.0/chain_madgraph.GetEntries();
	
	for(unsigned int i=0; i<chain_madgraph.GetEntries(); ++i)
	  {
	    if(i%10000==0)
	      cout<<"\r\033[K"<<floor((i/(chain_madgraph.GetEntries()*12.*3)  +
				       (3*(k-1)+ka)/(12.*3)
				       )*10000)*0.01
		  <<"% pre-processed (for prescale)..."
		  <<flush;
	  

	    chain_madgraph.GetEntry(i);
	    if(size<2) continue;

	    w*=rescale;
	    
	    int real_njet=0;
	    for(int n=0; n<size; ++n)
	      {
		if(t1[n] < 1e-5) continue;
		if(t2[n] < 1e-5) continue;
		if(pt[n] < min_pt) continue;
		if(m[n]<20) continue;
		real_njet++;
	      }
	    
	    int n_real=0;

	    //limit to two jets
	    for(int n=0; n<size && n_real<2; ++n)
	      {
		if(t1[n] < 1e-5) continue;
		if(t2[n] < 1e-5) continue;
		if(pt[n] < min_pt) continue;
		if(m[n]<20) continue;
		
		
		//double logpt=log(pt[n]/min_pt);
		//double mpt=-log10(m[n]/pt[n]);
		//double t21=t2[n]/t1[n];
		

		//see which prescale bin it is in
		if(n_real==2-1)
		  {
		    for(unsigned int il=0; il<3; il++)
		      {
			if(pt[n]> lhc_prescale_cut[il])
			  lhc_prescale[il].Add(w);
		      }
		  }

		n_real++;
	      }

	  }//loop over all entries	
      }//loop over all files 
  cout<<"\r\033[K100% pre-processed."<<endl<<flush;

  //now compute all the prescale factors
  for(unsigned int il=0; il<3; il++)
    {
      //the number correspond to cxn in millibarns
      //first compute event rate
      //assume instantenous luminosity = 10^{-34} cm^{-2}
      // = 10 nanobarn^{-1} = 1e7 milibarn^{-1}
      double event_rate= 1e7 * lhc_prescale[il].nevt();

      cout<<"event rate factor for PT bin "<<lhc_prescale_cut[il]
	  <<" is: "<<lhc_prescale[il].nevt()<<endl;

      //make sure event_rate is positive
      assert(event_rate>0);

      //we now want to resctrict it to 200 Hz
      lhc_prescale_factor[il]= lhc_Hz/ event_rate;

      if(lhc_prescale_factor[il]>1.0)
	lhc_prescale_factor[il]=1.0;

      //make sure it is positive
      cout<<"Prescale factor for PT bin "<<lhc_prescale_cut[il]
	  <<" is: "<<lhc_prescale_factor[il]<<endl;
      assert(lhc_prescale_factor[il]>0);
      
      
    }


}//end prescale
*/
  double total_cxn=0;

  for(unsigned int k=1; k<=12; k++)
    for(unsigned int ka=0; ka<3; ka++)
      {
	TChain chain_madgraph("Ntuple");
	stringstream sout;
	sout<<argv[1]<<alist[ka]<<"."<<k<<".root";
	chain_madgraph.Add(sout.str().c_str());
	
	//read the input file
	chain_madgraph.SetBranchAddress("size", &size);
	chain_madgraph.SetBranchAddress("w", &w);
	chain_madgraph.SetBranchAddress("pt", pt);
	chain_madgraph.SetBranchAddress("eta", eta);
	chain_madgraph.SetBranchAddress("phi", phi);
	chain_madgraph.SetBranchAddress("m", m);
	chain_madgraph.SetBranchAddress("tau1", t1);
	chain_madgraph.SetBranchAddress("tau2", t2);
	chain_madgraph.SetBranchAddress("tau3", t3);
	
	double rescale=1.0/chain_madgraph.GetEntries();
	
	for(unsigned int i=0; i<chain_madgraph.GetEntries(); ++i)
	  {
	    if(i%10000==0)
	      cout<<"\r\033[K"<<floor((i/(chain_madgraph.GetEntries()*12.*3)  +
				       (3*(k-1)+ka)/(12.*3)
				       )*10000)*0.01
		  <<"% trained..."
		  <<flush;
	  
	    
	    chain_madgraph.GetEntry(i);
	    if(size<2) continue;
	    
	    w*=rescale;

	    int real_njet=0;

	    //indices for the two leading jet
	    int n1,n2;
	    for(int n=0; n<size; ++n)
	      {
		if(t1[n] < 1e-5) continue;
		if(t2[n] < 1e-5) continue;
		if(t3[n] < 1e-5) continue;
		if(pt[n] < min_pt) continue;
		if(m[n]<20) continue;

		real_njet++;

		if(real_njet==1)
		  n1=n;
		if(real_njet==2)
		  {
		    n2=n;
		    break;
		  }
	      }
	    
	    if(real_njet<2) continue;


	    //first find out prescale factor
	    /*

	    for(unsigned int il=0; il<3; il++)
	      {
		//find prescale bin
		if(pt[n2]>= lhc_prescale_cut[il] && 
		   (il>=2 || pt[n2]< lhc_prescale_cut[il+1]))
		  {
		    double prescale_factor=lhc_prescale_factor[il];
		    //rescale the event
		    w*=prescale_factor;
		    break;
		  }
	      }
	    */

	    double logpt1=log(pt[n1]/min_pt);
	    double mpt1=-log10(m[n1]/pt[n1]);
	    double t32_1=t3[n1]/t2[n1];

	    double logpt2=log(pt[n2]/min_pt);
	    double mpt2=-log10(m[n2]/pt[n2]);
	    double t32_2=t3[n2]/t2[n2];


	    total_cxn+=w;		

	    if(!fast_mode)
	      data_m12->Fill(make_list<double>(m[n1], m[n2], t32_1, t32_2), w);
	    
	    //if one jet has mass < 140, fill template
	    if(m[n1]<mass_cut)
	      {
		
		if(!fast_mode)
		  {
		    dataset_mpt_njet[1]->Fill(make_list<double>(mpt2,logpt2),w);
		    dataset_tau32_njet[1]->Fill(make_list<double>(t32_2,logpt2),w);
		  }
		
		dataset_mpt_tau32_njet[0]->
		  Fill(make_list<double>(mpt2,t32_2,logpt2), w);


		//fill LHC statistics
		for(int is=0; is<lhc_size; ++is)
		  {
		    if(logpt2< is+1 && logpt2 > is)
		      lhc_stat[is]->Fill(mpt2, t32_2, w);
		  }

	      }


	    if(m[n2]<mass_cut)
	      {
		if(!fast_mode)
		  {
		    dataset_mpt_njet[0]->Fill(make_list<double>(mpt1,logpt1),
					      w);
		    dataset_tau32_njet[0]->Fill(make_list<double>(t32_1,logpt1),
						w);
		  }

		dataset_mpt_tau32_njet[0]->
		  Fill(make_list<double>(mpt1,t32_1,logpt1), w);

		//fill LHC statistics
		for(int is=0; is<lhc_size; ++is)
		  {
		    if(logpt1< is+1 && logpt1 > is)
		      lhc_stat[is]->Fill(mpt1, t32_1, w);
		  }


	      }

	    
	  }//loop over all entries
	
      }//loop over all files
  cout<<"\r\033[K100% trained."<<endl<<flush;


  //statistics plots
  TFile file_stat("output_stat.root", "RECREATE");
  
  cout<<"Statistics plots stored in output_stat.root"<<endl;
  for(int is=0; is<lhc_size; ++is)
    {
      //scale it up and get the weights
      for(int i=0; i<=lhc_stat[is]->GetNbinsX(); ++i)
      for(int j=0; j<=lhc_stat[is]->GetNbinsY(); ++j)
	{
	  double temp_scale=0;
	  if(lhc_stat[is]->GetBinContent(i,j) > 0)
	    temp_scale=lhc_factor*
	      lhc_stat[is]->GetBinError(i,j)*lhc_stat[is]->GetBinError(i,j)/
	      (lhc_stat[is]->GetBinContent(i,j));
	  lhc_stat[is]->SetBinContent(i, j, temp_scale);
	}

      lhc_stat[is]->Write();
    }

  file_stat.Close();

  //get PDF list
  vector<PDF> PDF_mpt_njet;
  vector<PDF> PDF_tau32_njet;
  vector<PDF> PDF_mpt_tau32_njet;
  
  //first turn-on Bias only in PDF computation
  DataSet::b_bias_only = true;
  
  //for the njet case
  for(int n=0; n<dataset_mpt_njet.size(); ++n)
    {
      if(!fast_mode)
	{
	  PDF_mpt_njet.push_back(dataset_mpt_njet[n]->ComputePDF());
	  PDF_tau32_njet.push_back(dataset_tau32_njet[n]->ComputePDF());      
	}
      //if PDF smearing is on, smear the pdf first
      if(b_smear_PDF)
	{
	  //if second jet, the PDF is also computed
	  if(n==1)
	    continue;

	  //a negative scale_factor means no scaling
	  double scale_factor=-1.0;

	  if(b_smear_lhc)
	    scale_factor=lhc_factor;
	  
	  DataSet dataset_new=
	    dataset_mpt_tau32_njet[n]->GenerateDataSet(scale_factor);
	  
	  PDF_mpt_tau32_njet.push_back(dataset_new.ComputePDF());
	  
	  //for first and second jet
	  if(n==0)
	      PDF_mpt_tau32_njet.push_back(dataset_new.ComputePDF());
	    
	}

      else
	PDF_mpt_tau32_njet.push_back(dataset_mpt_tau32_njet[n]->ComputePDF());
    }


  //free up memory
  for(unsigned int n=0; n<dataset_mpt_njet.size(); ++n)
    {
      if(!dataset_mpt_njet[n])
	{
	  delete dataset_mpt_njet[n];
	  dataset_mpt_njet[n]=NULL;
	}	  
    }

  for(unsigned int n=0; n<dataset_tau32_njet.size(); ++n)
    {
      if(!dataset_tau32_njet[n])
	{
	  delete dataset_tau32_njet[n];
	  dataset_tau32_njet[n]=NULL;
	}	  
    }

  for(unsigned int n=0; n<dataset_mpt_tau32_njet.size(); ++n)
    {
      if(!dataset_mpt_tau32_njet[n])
	{
	  delete dataset_mpt_tau32_njet[n];
	  dataset_mpt_tau32_njet[n]=NULL;
	}
    }


  //create generators for MC toy integration
  vector<PDF::MC_Generator> MC_mpt_tau32_njet;
  
  for(int n=0; n< PDF_mpt_tau32_njet.size(); ++n)
    {
      MC_mpt_tau32_njet.push_back(PDF_mpt_tau32_njet[n].Generator
				  (make_list<PDF::MC_Flag>
				   (PDF::OUTPUT, PDF::OUTPUT, PDF::INPUT)
				   ));
    }


  //declare 2D distributions to study
  DataSet t32_j12(make_list<double>(0, 0),
		  make_list<double>(1.0, 1.0),
		  make_list<double>(0.02, 0.02));
  
  DataSet mpt_j12(make_list<double>(0, 0),
		  make_list<double>(1.5, 1.5),
		  make_list<double>(0.01, 0.01));

  DataSet t32_j12_MC(make_list<double>(0, 0),
		     make_list<double>(1.0, 1.0),
		     make_list<double>(0.02, 0.02));
  
  DataSet mpt_j12_MC(make_list<double>(0, 0),
		     make_list<double>(1.5, 1.5),
		     make_list<double>(0.01, 0.01));

  TH2D t32_12("t32_12", "2nd jet #tau_{32} vs. first jet #tau_{32}", 
	   200, 0, 1, 200, 0, 1);

  TH2D t32_12_MC("t32_12_MC", "2nd jet #tau_{32} vs. first jet #tau_{32} (MC)", 
		 200, 0, 1, 200, 0, 1);


  TH2D m12_3("m12_3", "2nd jet mass vs. first jet mass (0.3)", 
	   170, 300, 2000, 170, 300, 2000);

  TProfile2D m12_smooth_3("m12_smooth_3", "2nd jet mass vs. first jet mass (0.3 smooth)", 
	   170, 300, 2000, 170, 300, 2000);

  TH2D m12_MC_3("m12_MC_3", "2nd jet mass vs. first jet mass (MC) (0.3)", 
	      170, 300, 2000, 170, 300, 2000);

  TH2D m12_4("m12_4", "2nd jet mass vs. first jet mass (0.4)", 
	   170, 300, 2000, 170, 300, 2000);

  TProfile2D m12_smooth_4("m12_smooth_4", "2nd jet mass vs. first jet mass (0.4 smooth)", 
	   170, 300, 2000, 170, 300, 2000);

  TH2D m12_MC_4("m12_MC_4", "2nd jet mass vs. first jet mass (MC) (0.4)", 
	      170, 300, 2000, 170, 300, 2000);

  TH2D m12_5("m12_5", "2nd jet mass vs. first jet mass (0.5)", 
	   170, 300, 2000, 170, 300, 2000);

  TProfile2D m12_smooth_5("m12_smooth_5", "2nd jet mass vs. first jet mass (0.5 smooth)", 
	   170, 300, 2000, 170, 300, 2000);

  TH2D m12_MC_5("m12_MC_5", "2nd jet mass vs. first jet mass (MC) (0.5)", 
	      170, 300, 2000, 170, 300, 2000);

  TH2D m12_6("m12_6", "2nd jet mass vs. first jet mass (0.6)", 
	   170, 300, 2000, 170, 300, 2000);

  TProfile2D m12_smooth_6("m12_smooth_6", "2nd jet mass vs. first jet mass (0.6 smooth)", 
	   170, 300, 2000, 170, 300, 2000);

  TH2D m12_MC_6("m12_MC_6", "2nd jet mass vs. first jet mass (MC) (0.6)", 
	      170, 300, 2000, 170, 300, 2000);



  if(!fast_mode)
    {

  cout<<"computing smoothed 2D histograms"<<endl;
  PDF m12_pdf=data_m12->ComputePDF();
  delete data_m12;

  PDF::MC_Generator m12_gen=m12_pdf.Generator(
   make_list<PDF::MC_Flag>(PDF::OUTPUT, PDF::OUTPUT, 
			   PDF::OUTPUT, PDF::OUTPUT)
					      );
  
  //get the full integration for 
  vector<PDF::MC_Event> m12_events=
    m12_gen.GenerateFull(vector<double>());


  for(unsigned int i=0; i<m12_events.size(); ++i)
    {
      double m1,m2,t32_1,t32_2;
      m1=m12_events[i].value[0];
      m2=m12_events[i].value[1];
      t32_1=m12_events[i].value[2];
      t32_2=m12_events[i].value[3];
      double ww=m12_events[i].weight * total_cxn;

      if(m1 < mass_cut || m2 < mass_cut) continue;

      for(unsigned int n=0; n<t32_cut.size(); n++)
	{
	  
	  if(t32_1 < 0.3 && 
	     t32_2 < 0.3 )
	    m12_smooth_3.Fill(m1, m2, ww);
	  
	  if(t32_1 < 0.4 && 
	     t32_2 < 0.4 )
	    m12_smooth_4.Fill(m1, m2, ww);
		
	  if(t32_1 < 0.5 && 
	     t32_2 < 0.5 )
	    m12_smooth_5.Fill(m1, m2, ww);
	  
	  if(t32_1 < 0.6 && 
	     t32_2 < 0.6 )
	    m12_smooth_6.Fill(m1, m2, ww);
	  
	}
    }

    }

  cout<<"starting toy MC integration..."<<endl;
  
  //declare event counter
  double n_evt=0;

  for(unsigned int k=1; k<=12; k++)
    for(unsigned int ka=1; ka<3; ka++)
      {
	TChain chain_madgraph("Ntuple");
	stringstream sout;
	sout<<argv[1]<<alist[ka]<<"."<<k<<".root";
	chain_madgraph.Add(sout.str().c_str());
	
	//read the input file
	chain_madgraph.SetBranchAddress("size", &size);
	chain_madgraph.SetBranchAddress("w", &w);
	chain_madgraph.SetBranchAddress("pt", pt);
	chain_madgraph.SetBranchAddress("eta", eta);
	chain_madgraph.SetBranchAddress("phi", phi);
	chain_madgraph.SetBranchAddress("m", m);
	chain_madgraph.SetBranchAddress("tau1", t1);
	chain_madgraph.SetBranchAddress("tau2", t2);
	chain_madgraph.SetBranchAddress("tau3", t3);
	
	double rescale=1.0/chain_madgraph.GetEntries();
	

	int init_entry=0;
	int final_entry=chain_madgraph.GetEntries();
	

	//if dividing up the jobs
	if(b_divide)
	  {
	    int div_size=(chain_madgraph.GetEntries()/n_div)+1;
	    init_entry=div_size*n_part;
	    final_entry=init_entry+div_size;

	    //make sure the index is less than the maximum
	    if(final_entry > chain_madgraph.GetEntries())
	      final_entry = chain_madgraph.GetEntries();
	  }

	

	for(unsigned int i=init_entry; i<final_entry; ++i)
	  {
	    if(i%10000==0)
	      cout<<"\r\033[K"<<floor((i/(chain_madgraph.GetEntries()*12.*3)  +
				       (3*(k-1)+ka)/(12.*3)
				       )*10000)*0.01
		  <<"% processed..."
		  <<flush;

	    chain_madgraph.GetEntry(i);
	    if(size<2) continue;

	    w*=rescale;
	    
	    //if doing lhc scaling, rescale the weight
	    if(b_smear_lhc)
	      w*=lhc_factor;

	    int real_njet=0;
	    int n1,n2;

	    for(int n=0; n<size; ++n)
	      {
		if(t1[n] < 1e-5) continue;
		if(t2[n] < 1e-5) continue;
		if(t3[n] < 1e-5) continue;

		//large pt cut
		if(pt[n] < 350.) continue;
		if(m[n]<20) continue;
		
		real_njet++;

		if(real_njet==1)
		  n1=n;
		if(real_njet==2)
		  {
		    n2=n;
		    break;
		  }
	      }
	    
	    //demand at least 2 jets
	    if(real_njet < 2) continue;
	    
	    double logpt1=log(pt[n1]/min_pt);
	    double mpt1=-log10(m[n1]/pt[n1]);
	    double t32_1=t3[n1]/t2[n1];

	    double logpt2=log(pt[n2]/min_pt);
	    double mpt2=-log10(m[n2]/pt[n2]);
	    double t32_2=t3[n2]/t2[n2];

	    //now populate the DataSet
	    t32_j12.Fill(make_list<double>(t32_1, t32_2), w*rescale);
	    mpt_j12.Fill(make_list<double>(mpt1, mpt2), w*rescale);


	    //fill the counters
	    for(unsigned int n=0; n<t32_cut.size(); n++)
	      {
		nevt_ABCD[n].Add(m[n1], m[n2],
				 t32_1, t32_2, w);
		
		if( t32_1 < t32_cut[n] &&
		    t32_2 < t32_cut[n] &&
		    m[n1] > mj_cut[n] &&
		    m[n2] > mj_cut[n] )		 
		  nevt_real[n].Add(w);		
	      }

	    
	    if(m[n1] > mass_cut && m[n2] > mass_cut)
	      {
		if(t32_1 < 0.3 && 
		   t32_2 < 0.3 )
		  m12_3.Fill(m[n1], m[n2], w);
		
		if(t32_1 < 0.4 && 
		   t32_2 < 0.4 )
		  {
		    m12_4.Fill(m[n1], m[n2], w);
		    t32_12.Fill(t32_1, t32_2, w);
		  }
		
		if(t32_1 < 0.5 && 
		   t32_2 < 0.5 )
		    m12_5.Fill(m[n1], m[n2],w);
		
		if(t32_1 < 0.6 && 
		   t32_2 < 0.6 )
		  m12_6.Fill(m[n1], m[n2],w);
	      }
	    //we have an event now!
	    n_evt+=w;

	    //now the first three jets, we have templates
	    //so create a vector for them
	    vector<vector<double> > MC_m;
	    vector<vector<double> > MC_mpt;
	    vector<vector<double> > MC_tau32;
	    vector<vector<double> > MC_weight;
	    vector<vector<double> > MC_weight_bias;
	    vector<vector<double> > MC_err;
	    
	    int nGenerate=500;

	    //now populate the MC data
	    vector<PDF::MC_Event> MC_mpt_tau32_1=
	      MC_mpt_tau32_njet[0].GenerateFast
	      (make_list<double>(log(pt[n1]/min_pt)), nGenerate);

	    vector<PDF::MC_Event> MC_mpt_tau32_2=
	      MC_mpt_tau32_njet[0].GenerateFast
	      (make_list<double>(log(pt[n2]/min_pt)), nGenerate);


	    //two jets, hence two vectors
	    MC_m.push_back(vector<double>(MC_mpt_tau32_1.size()));
	    MC_mpt.push_back(vector<double>(MC_mpt_tau32_1.size()));
	    MC_tau32.push_back(vector<double>(MC_mpt_tau32_1.size()));
	    MC_weight.push_back(vector<double>(MC_mpt_tau32_1.size()));
	    MC_weight_bias.push_back(vector<double>(MC_mpt_tau32_1.size()));
	    MC_err.push_back(vector<double>(MC_mpt_tau32_1.size()));

	    MC_m.push_back(vector<double>(MC_mpt_tau32_1.size()));
	    MC_mpt.push_back(vector<double>(MC_mpt_tau32_1.size()));
	    MC_tau32.push_back(vector<double>(MC_mpt_tau32_1.size()));
	    MC_weight.push_back(vector<double>(MC_mpt_tau32_1.size()));
	    MC_weight_bias.push_back(vector<double>(MC_mpt_tau32_1.size()));
	    MC_err.push_back(vector<double>(MC_mpt_tau32_1.size()));
	    
	    for(unsigned int nMC=0; nMC<MC_mpt_tau32_1.size(); ++nMC)
	      {
		double t_mass1=exp10(-MC_mpt_tau32_1[nMC].value[0])*pt[n1];	
		double t_mass2=exp10(-MC_mpt_tau32_2[nMC].value[0])*pt[n2];
		
		double t_tau32_1=MC_mpt_tau32_1[nMC].value[1];
		double t_tau32_2=MC_mpt_tau32_2[nMC].value[1];

		double t_mpt_1=MC_mpt_tau32_1[nMC].value[0];
		double t_mpt_2=MC_mpt_tau32_2[nMC].value[0];
		
		MC_m[0].push_back(t_mass1);
		MC_m[1].push_back(t_mass2);

		MC_tau32[0].push_back(t_tau32_1);
		MC_tau32[1].push_back(t_tau32_2);

		MC_mpt[0].push_back(t_mpt_1);
		MC_mpt[1].push_back(t_mpt_2);

		MC_weight[0].push_back(MC_mpt_tau32_1[nMC].weight);
		MC_weight[1].push_back(MC_mpt_tau32_2[nMC].weight);

		MC_weight_bias[0].push_back
		  (MC_mpt_tau32_1[nMC].weight-MC_mpt_tau32_1[nMC].err);

		MC_weight_bias[1].push_back
		  (MC_mpt_tau32_2[nMC].weight-MC_mpt_tau32_2[nMC].err);

		MC_err[0].push_back(MC_mpt_tau32_1[nMC].err);
		MC_err[1].push_back(MC_mpt_tau32_2[nMC].err);
	      }
	    
	    //first reweight the events to make the sum zero
	    double weight12=0;
	    double weight12_bias=0;

	    for(unsigned int n=0; n<MC_weight[0].size(); ++n)
	      {
		weight12 += MC_weight[0][n]*MC_weight[1][n];
		weight12_bias += MC_weight_bias[0][n]*MC_weight_bias[1][n];
	      }

	    assert(weight12!=0);
	    assert(weight12_bias!=0);

	    weight12=1.0/weight12;
	    weight12_bias=1.0/weight12_bias;
	    
	    //now we have a bunch of toy MC events
	    //populate the DataSet using MC events
	    //and get the fraction of events that pass
	    //first we create temporary vector to store error
	    vector<double> nevt_MC_temp(nevt_MC.size());
	    vector<double> nevt_err_temp(nevt_MC.size());

	    for(unsigned int n=0; n<MC_weight[0].size(); ++n)
	      {
		//now populate the DataSet
		t32_j12_MC.Fill(make_list<double>
				(MC_tau32[0][n], MC_tau32[1][n]),
				w*weight12*
				MC_weight[0][n]*MC_weight[1][n]);
		
		mpt_j12_MC.Fill(make_list<double>
				(MC_mpt[0][n], MC_mpt[1][n]), 
				w*weight12*
				MC_weight[0][n]*MC_weight[1][n]);
		
		double weight=MC_weight[0][n]*MC_weight[1][n]*weight12;

		double err= MC_weight_bias[0][n]*MC_weight_bias[1][n]*
		  weight12_bias;



		if(MC_m[0][n] > mass_cut && MC_m[1][n] > mass_cut)
		  {
		    
		    if(MC_tau32[0][n] < 0.3 && 
		       MC_tau32[1][n] < 0.3 )
			m12_MC_3.Fill(MC_m[0][n], MC_m[1][n], err*w);

		    if(MC_tau32[0][n] < 0.4 && 
		       MC_tau32[1][n] < 0.4 )
		      {
			t32_12_MC.Fill(MC_tau32[0][n], MC_tau32[1][n], err*w);
			m12_MC_4.Fill(MC_m[0][n], MC_m[1][n], err*w);
		      }

		    if(MC_tau32[0][n] < 0.5 && 
		       MC_tau32[1][n] < 0.5 )
		      m12_MC_5.Fill(MC_m[0][n], MC_m[1][n], err*w);
		    
		    if(MC_tau32[0][n] < 0.6 && 
		       MC_tau32[1][n] < 0.6 )
		      m12_MC_6.Fill(MC_m[0][n], MC_m[1][n], err*w);
		  }

		for(unsigned int nn=0; nn<t32_cut.size(); nn++)
		  {
		    if(MC_m[0][n] > mj_cut[nn] && MC_m[1][n] > mj_cut[nn] && 
		       MC_tau32[0][n] < t32_cut[nn] && 
		       MC_tau32[1][n] < t32_cut[nn])
		      {
			//use the corrected weight
			nevt_MC_temp[nn] += weight;
			nevt_err_temp[nn] += err;
		      }
		  }

	      } // loop over all MC events
	    
	    //increment count for data-driven
	    for(unsigned int n=0; n<t32_cut.size(); n++)
	      {
		nevt_MC[n].Add(w*nevt_MC_temp[n],
			       w*(nevt_MC_temp[n]-nevt_err_temp[n]));
	      }


	  }//loop all entries in file
      }//loop over all madgraph files
  cout<<"\r\033[Kfinished all processing"<<endl<<flush;
  

  cout<<"Writing graphs..."<<endl;

  //store some of the templates
  TString outfilename(argv[2]);
  if(b_divide)
    {
      outfilename += "_job";
      outfilename += n_part;      
    }
  outfilename += ".root";
  TFile out(outfilename, "RECREATE");

  plot2D(t32_j12, "t32_j12").Write();
  plot2D(t32_j12_MC, "t32_j12_MC").Write();

  plot2D(mpt_j12, "mpt_j12").Write();
  plot2D(mpt_j12_MC, "mpt_j12_MC").Write();

  out.mkdir("templates")->cd();
  map<string, Vec> graph_list;


  if(!fast_mode)
  for(int i=60; i<=3000; i+=10)
    {
      stringstream sout;
      
      //do the njet case
      for(int n=0; n<2; ++n)
	{
	  sout.str("");
	  sout<<"t32_"<<n<<"jet_"<<i;
	  graph_list[sout.str()].fill_data(PDF_tau32_njet[n], i);
	  graph_list[sout.str()].get_graph(sout.str()).Write();

	  sout.str("");
	  sout<<"mpt_"<<n<<"jet_"<<i;
	  graph_list[sout.str()].fill_data(PDF_mpt_njet[n], i);
	  graph_list[sout.str()].get_graph(sout.str(), true).Write();	  
	}

    }
  
  m12_3.Write();
  m12_MC_3.Write();
  m12_smooth_3.Write();
  m12_4.Write();
  m12_MC_4.Write();
  m12_smooth_4.Write();
  m12_5.Write();
  m12_MC_5.Write();
  m12_smooth_5.Write();
  m12_6.Write();
  m12_MC_6.Write();
  m12_smooth_6.Write();

  t32_12.Write();
  t32_12_MC.Write();  

  out.Close();

  TString cuttablename("cut_table");

  if(b_divide)
    {
      cuttablename += "_job";
      cuttablename += n_part;    
    }

  cuttablename+=".txt";
    
  ofstream fout(cuttablename);

  fout<<"t32_cut, mj_cut, nevt_real, nevt_real_err, nevt_MC, nevt_MC_bias, stat_err, alpha, alpha_err, n_kinematic"<<endl;

  for(unsigned int n=0; n<t32_cut.size(); n++)
    {
      fout<<t32_cut[n]<<", "
	  <<mj_cut[n]<<", "
	  <<nevt_real[n].nevt()<<", "<<nevt_real[n].err()<<", "
	  <<nevt_MC[n].nevt()<<", "
	  <<nevt_MC[n].bias()<<", "
	  <<nevt_MC[n].var()<<", "
	  <<nevt_ABCD[n].factor()<<", "
	  <<nevt_ABCD[n].factor_err()
	  <<endl;
    }

  fout<<endl<<endl<<endl<<endl;
  fout<<"output corrected results"<<endl;
  fout<<"t32_cut,mj_cut,n_corrected"<<endl;
  for(unsigned int n=0; n<t32_cut.size(); n++)
    {
      fout<<t32_cut[n]<<", "
	  <<mj_cut[n]<<", "
	  <<nevt_MC[n].nevt_corr()
	  <<endl;
    }


  fout.close();

  cout<<"Kinematic sample: "<<n_evt<<endl;
  cout<<"Program ends."<<endl;
  return 1;
}
  
